home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / raytrace / renaisnc / rcnstrct.lha / Reconstruct / Reconstruct.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-01-04  |  24.2 KB  |  928 lines
  1. /**
  2.  **  author    dan O'Donnell, James Painter
  3.  **  purpose   Reconstruct the image function from stochastic sample file 
  4.  **
  5.  **/
  6.  
  7. #include <stdio.h>
  8. #include <math.h>
  9. #include <assert.h>
  10. #include "Reconstruct.h"
  11. #include "wff.h" 
  12.  
  13. /*
  14. ** Test to see if two rectangles overlap.  Return TRUE if the do, FALSE
  15. ** otherwise.
  16. */
  17. #define RectOverlap(Rect1,Rect2)  \
  18.  ( ! ( (Rect1)->xH < (Rect2)->xL || (Rect2)->xH < (Rect1)->xL ||     \
  19.        (Rect1)->yH < (Rect2)->yL || (Rect2)->yH < (Rect1)->yL    ))  
  20.  
  21.  
  22. /*------------------------------------
  23.   Check to see if the point X,Y is in 
  24.   the rectangle Rect.
  25.  ------------------------------------*/
  26. #define PointInRectangle(Xs,Ys,Rect)                     \
  27.     ( ((Xs) >= (Rect)->xL)  &&  ((Xs) < (Rect)->xH)  &&  \
  28.       ((Ys) >= (Rect)->yL)  &&  ((Ys) < (Rect)->yH) )
  29.  
  30.  
  31.  
  32. /*  Number of nodes/samples to allocate at a time */
  33. #define GLOB_SIZE 1024
  34.  
  35.  
  36.  
  37. /*----------------------------------------
  38.   Subdivision level
  39.   --------------------------------------*/
  40. static
  41. int level;
  42.  
  43. int AdaptiveSplits = FALSE;
  44. int Animating = FALSE;
  45.  
  46. HierarchicalRegion Root;
  47. RectType RootBoundary;
  48.  
  49. extern double atof();
  50. extern  HierarchicalRegion *PropagateSample();
  51. extern  FloatColor FilterFinalize();
  52. extern  char *mallocNd();
  53.  
  54. int window;            /* TRUE if sinc function is windowed */
  55.  
  56. static  FILE *NumRaysImage = NULL;
  57. int Xsize = 128, Ysize = 128;
  58. int x_start=0, y_start=0, x_stop=127, y_stop=127;
  59. double PixelSizeX, PixelSizeY;
  60. double MaxInBandValue = 1.0; /* The maximum value expected in a color band */
  61.  
  62. double FilterSupport = 3.0;
  63. static  double Aspect = (6.0/5.3), FilterRadius = 1.25;
  64. double B = .3333333333, C = .33333333333; /* "Mitchell" cubic parameters */
  65. double CutOff = -1.0;
  66. FilterFunctionType FilterFunction = Gaussian;
  67.  
  68. int BitsPerBand = 8;        /* Output bits per band */
  69. long MaxOutBandValue = 255; /* The maximum value to output */
  70.  
  71. FILE *image = NULL;
  72. char *filename;
  73.  
  74. typedef struct {
  75.   FloatColor clr;
  76.   float weight;
  77.   unsigned short samplesInPixel;
  78. } Pixel;
  79.  
  80. Pixel **image_data;
  81.  
  82. /* --------------------------------------------------------------------------
  83.    Split a node boundary rectangle into the rectangle for the Left 
  84.    (if LeftOrRight == TRUE) or Right (if LeftOrRight == FALSE) subnode.
  85.    ---------------------------------------------------------------------------*/
  86. SplitBoundary(LeftOrRight, splitDirection, splitValue,
  87.           OldBoundary, NewBoundary)
  88.      int LeftOrRight;
  89.      double splitValue;
  90.      RectType *OldBoundary, *NewBoundary;
  91. {
  92.   double Xs, Ys;
  93.   
  94.   *NewBoundary = *OldBoundary;
  95.   if ( splitDirection )  /* Odd levels split in X */
  96.     {  /* Split in X */
  97.       Xs = splitValue;
  98.       if (LeftOrRight) {    /* Left Child */
  99.     NewBoundary->xH = Xs;
  100.       } else {            /* Right child */
  101.     NewBoundary->xL = Xs; 
  102.       }
  103.     }
  104.   else
  105.     {  /* Split in Y */
  106.       Ys = splitValue;
  107.       if (LeftOrRight) {    /* Left Child  */
  108.     NewBoundary->yH = Ys;
  109.       } else {            /* Right Child */
  110.     NewBoundary->yL = Ys; 
  111.       }
  112.     }
  113. }
  114.  
  115.  
  116.  
  117. /*-----------------------
  118.   allocate a node
  119.   ----------------------*/
  120.  
  121. static
  122.   HierarchicalRegion *
  123.   CreateNode()
  124. {
  125.   char *calloc();
  126.   static HierarchicalRegion *FreeNodes;
  127.   static nFreeNodes = 0;
  128.  
  129.  
  130.   if (nFreeNodes <= 0) {
  131.     FreeNodes = (HierarchicalRegion *) 
  132.                   calloc( GLOB_SIZE, sizeof(HierarchicalRegion));
  133.     if(FreeNodes == NULL)
  134.       {
  135.     fprintf(stderr,"Fatal error - Can't allocate nodes\n");
  136.     abort();
  137.       }
  138.     nFreeNodes = GLOB_SIZE;
  139.   }
  140.   return &(FreeNodes[--nFreeNodes]);
  141. }
  142.  
  143.  
  144. static SampleData *SampleFreeList = NULL;
  145.  
  146. FreeSample( Sample )
  147. SampleData *Sample;
  148. {
  149.   /* Link the sample onto the free list */
  150.   *( (SampleData **) Sample) = SampleFreeList;
  151.   SampleFreeList = Sample;
  152. }
  153.  
  154. SampleData *CreateSample()
  155. {
  156.   char *calloc();
  157.   static SampleData *FreeSamples, *Sample;
  158.   static int nFreeSamples = 0;
  159.  
  160.   /* Get one off the free list if possible */
  161.   if (SampleFreeList != NULL) {
  162.     Sample = SampleFreeList;
  163.     SampleFreeList = *( (SampleData **) SampleFreeList );
  164.     bzero ( Sample, sizeof(SampleData) );
  165.     return Sample;
  166.   }
  167.  
  168.   if (nFreeSamples <= 0) {
  169.     FreeSamples = (SampleData *) calloc( GLOB_SIZE, sizeof(SampleData));
  170.     if(FreeSamples == NULL)
  171.       {
  172.     fprintf(stderr,"Fatal error - Can't allocate sample data\n");
  173.     abort();
  174.       }
  175.     nFreeSamples = GLOB_SIZE;
  176.   }
  177.   return & (FreeSamples[--nFreeSamples]);
  178. }
  179.  
  180. /* Split a node into two subnodes */
  181. SplitNode(Node)
  182.      HierarchicalRegion *Node;
  183. {
  184.   Node->Left = CreateNode();
  185.   Node->Right = CreateNode();
  186. }
  187.  
  188.  
  189. /*-----------------------------
  190.   propagate the current node's
  191.   sample to one of the subnodes
  192.   return the EMPTY SUBNODE
  193.   ---------------------------*/
  194. HierarchicalRegion *PropagateSample(Node, Boundary, NewSample)
  195. HierarchicalRegion *Node;
  196. RectType *Boundary;
  197. SampleData *NewSample;
  198. {
  199.   RectType LeftBoundary;
  200.   double dx, dy;
  201.  
  202.   if (AdaptiveSplits) {
  203.  
  204.     /** Choose a split direction which separates the new sample
  205.      ** from the old and maintains an aspect ratio as close to square
  206.      ** as possible for the new sub cells.
  207.      */
  208.     double xsplit, ysplit, x1, x2, y1, y2, xratio, xratio2, yratio, yratio2;
  209.     double width, height;
  210.     
  211.  
  212.     /* Split at the mid point between samples */
  213.     xsplit = 0.5*(Node->Sample->Xs + NewSample->Xs);
  214.     ysplit = 0.5*(Node->Sample->Ys + NewSample->Ys);
  215.  
  216.     /* These are the side lengths of the sub cells */
  217.     x1 = xsplit - Boundary->xL;
  218.     x2 = Boundary->xH - xsplit;
  219.     y1 = ysplit - Boundary->yL;
  220.     y2 = Boundary->yH - ysplit;
  221.  
  222.     /* These are the side lengths of the whole cell */
  223.     width  = (Boundary->xH - Boundary->xL);
  224.     height = (Boundary->yH - Boundary->yL);
  225.  
  226.     /* These are the aspect ratios if we split in x */
  227.     if (x1 > height)
  228.       xratio =  (height < 1.0E-20) ? 1.0E20 : x1 / height;
  229.     else
  230.       xratio =  (x1     < 1.0E-20) ? 1.0E20 : height / x1;
  231.     if (x2 > height)
  232.       xratio2 =  (height < 1.0E-20) ? 1.0E20 : x2 / height;
  233.     else
  234.       xratio2 =  (x2     < 1.0E-20) ? 1.0E20 : height / x2;
  235.  
  236.     
  237.     /* Use the worst of the two for comparison */
  238.     if (xratio2 > xratio)  xratio = xratio2;
  239.  
  240.     /* These are the aspect ratios if we split in y */
  241.     if (y1 > width)
  242.       yratio =  (width < 1.0E-20) ? 1.0E20 : y1 / width;
  243.     else
  244.       yratio =  (y1     < 1.0E-20) ? 1.0E20 : width / y1;
  245.     if (y2 > width)
  246.       yratio2 =  (width < 1.0E-20) ? 1.0E20 : y2 / width;
  247.     else
  248.       yratio2 =  (y2     < 1.0E-20) ? 1.0E20 : width / y2;
  249.  
  250.  
  251.     /* Use the worst of the two for comparison */
  252.     if (yratio2 > yratio)  yratio = yratio2;
  253.  
  254.     if (xratio < yratio) {  /* X split is better */
  255.       Node->splitDirection = 1;    /* Split in X */
  256.       Node->splitValue = 0.5*(Node->Sample->Xs + NewSample->Xs);
  257.     } else  {
  258.       Node->splitDirection = 0;    /* Split in Y */
  259.       Node->splitValue = 0.5*(Node->Sample->Ys + NewSample->Ys);
  260.     }
  261.   } else {
  262.     Node->splitDirection = (level & 1);  /* Odd levels split in X */
  263.     if (Node->splitDirection)
  264.       Node->splitValue = 0.5*(Boundary->xL + Boundary->xH);
  265.     else
  266.       Node->splitValue = 0.5*(Boundary->yL + Boundary->yH);
  267.   }
  268.     
  269.   SplitBoundary( TRUE, Node->splitDirection, Node->splitValue,
  270.           Boundary, &LeftBoundary );
  271.   
  272.   if( (Node->Sample->Xs >=   LeftBoundary.xL)  &&
  273.       (Node->Sample->Xs <    LeftBoundary.xH)  &&
  274.       (Node->Sample->Ys >=   LeftBoundary.yL)  &&
  275.       (Node->Sample->Ys <    LeftBoundary.yH) )
  276.     {
  277.       Node->Left->Sample    = Node->Sample;
  278.       return Node->Right; 
  279.     }
  280.   else
  281.     {
  282.       Node->Right->Sample    = Node->Sample;
  283.       return Node->Left;
  284.     } 
  285. }
  286.  
  287.  
  288. SampleData *ReadSample(file)
  289. FILE *file;
  290. {  
  291.   SampleData *Sample = CreateSample();
  292.   long sibling;            /* Index of sibling node */
  293.  
  294. #ifdef BINARY
  295.   {
  296.     long array[7], *p = array;
  297.     if (fread( array, sizeof(long), 7, file ) != 7) return 0;
  298.     sibling    = *((long *) p++);
  299.     Sample->Xs = *((float *)p++);
  300.     Sample->Ys = *((float *)p++);
  301.     Sample->Ave.r = *((float *)p++);
  302.     Sample->Ave.g = *((float *)p++);
  303.     Sample->Ave.b = *((float *)p++);
  304.     Sample->Ave.a = *((float *)p);
  305.   }
  306. #else
  307.   if (fscanf( file, "%d %f %f %f %f %f %f", 
  308.            &sibling, &Sample->Xs, &Sample->Ys, 
  309.             &Sample->Ave.r, &Sample->Ave.g, &Sample->Ave.b, &Sample->Ave.a )
  310.       
  311.       != 7) return 0;
  312. #endif
  313.   return Sample;
  314. }
  315.  
  316.  
  317.  
  318.  
  319. /*---------------------------------
  320.   Insert the new node in the tree.
  321.  ----------------------------------*/
  322. InsertNode(NewSample, Root, RootBoundary, InterestingBoundary)
  323.      SampleData *NewSample;
  324.      HierarchicalRegion *Root;
  325.      RectType *RootBoundary, *InterestingBoundary;
  326. {
  327.   HierarchicalRegion *Node = Root, *EmptyNode;
  328.   RectType Boundarys[2], *Boundary, *NewBoundary;
  329.   int old, new;
  330.   
  331.  
  332.   Boundarys[0] = *RootBoundary;
  333.   old = 0; new = 1;
  334.   Boundary = Boundarys+0;
  335.   NewBoundary = Boundarys+1;
  336.  
  337.   /* Search down to the leaf level for the cell this sample belongs in */
  338.   level = 0;
  339.   while (!IsLeaf(Node))
  340.     {
  341.       if (!RectOverlap( InterestingBoundary, Boundary )) {
  342.     /* Not in an interesting part of the image, ignore it */
  343.     FreeSample( NewSample );
  344.     return;
  345.       }
  346.       SplitBoundary( TRUE, Node->splitDirection, Node->splitValue,
  347.             Boundary, NewBoundary );
  348.       if (PointInRectangle( NewSample->Xs, NewSample->Ys, NewBoundary )) {
  349.     Node = Node->Left;
  350.       } else {
  351.     SplitBoundary( FALSE, Node->splitDirection, Node->splitValue,
  352.               Boundary, NewBoundary );
  353.     Node = Node->Right;
  354.       }
  355.       old = (old +1) %2;
  356.       new = (new +1) %2;
  357.       Boundary    = Boundarys+old;
  358.       NewBoundary = Boundarys+new;
  359.       level++;
  360.     }
  361.   
  362.   /* Split the Node */
  363.   SplitNode( Node );
  364.  
  365.   /* Fill in the two new cells */
  366.   EmptyNode = PropagateSample(Node, Boundary, NewSample);
  367.   EmptyNode->Sample = NewSample;
  368.  
  369.   if (Animating)
  370.     AnimateShadeCell( Boundary, NewSample );
  371. }
  372.  
  373.  
  374. /*-----------------------
  375.   Read the samples
  376.   ----------------------*/
  377. ReadSampleFile( samplefile, InterestingBoundary )
  378. FILE *samplefile;
  379. RectType *InterestingBoundary;    /*  Region of interest  */
  380. {
  381.   /* read the boundary Rectangle */
  382. #ifdef BINARY
  383.   {
  384.     float array[4], *p=array;
  385.     fread( array, sizeof(float), 4, samplefile );
  386.     RootBoundary.xL = *p++;
  387.     RootBoundary.yL = *p++;
  388.     RootBoundary.xH = *p++;
  389.     RootBoundary.yH = *p;
  390.   }
  391. #else
  392.   fscanf( samplefile, "%f %f %f %f", 
  393.      &RootBoundary.xL, &RootBoundary.yL, 
  394.      &RootBoundary.xH, &RootBoundary.yH );
  395. #endif
  396.  
  397.   Root.Left = Root.Right = (HierarchicalRegion *) NULL;
  398.  
  399.   /* Read the Root Node */
  400.   Root.Sample = ReadSample( samplefile );
  401.  
  402.   while (!feof(samplefile)) {
  403.     SampleData *Sample;
  404.  
  405.     if (Animating)
  406.       AnimateCheckInput();
  407.  
  408.     /* Read the next node */
  409.     Sample = ReadSample( samplefile );
  410.     if (Sample == NULL)
  411.       break;
  412.  
  413.     /* Insert this node in the tree */
  414.     InsertNode( Sample, &Root, &RootBoundary, InterestingBoundary );
  415.   }
  416. }
  417.  
  418.  
  419.  
  420.  
  421. /*
  422. **   Traverse the sample tree and add the contribution from each sample
  423. **   into all affected pixels.
  424. */
  425. static
  426. Filter( Root, RootBoundary )
  427. HierarchicalRegion *Root;
  428. RectType RootBoundary;
  429. {
  430.   register HierarchicalRegion *Node;
  431.   RectType LeftBoundary, RightBoundary, Boundary, Overlap;
  432.   int Left, Right;
  433.   float Xs, Ys;
  434.  
  435.   /* This stack is used to avoid the expense of recursion */
  436.   HierarchicalRegion *NodeStack[MAX_DEPTH];
  437.   RectType BoundaryStack[MAX_DEPTH];
  438.   int levelStack[MAX_DEPTH];
  439.   int nStack;
  440.   int r, c, rmn, rmx, cmn, cmx;
  441.   REAL x, y, yL, yH, xL, xH, Weight;
  442.   Pixel *pixel;
  443.   FloatColor *clr;
  444.   extern double Gauss(), WindowedSincFunction(), CdfLookup();
  445.  
  446.   switch (FilterFunction) {
  447.   case Gaussian:
  448.     SetFilterSupport( FilterSupport );
  449.     BuildCdf( Gauss );
  450.     break;
  451.   case Sinc:
  452.     FilterSupport = FilterRadius;
  453.     SetFilterSupport( FilterSupport );
  454.     window = 0;
  455.     BuildCdf( WindowedSincFunction );
  456.     break;
  457.   case WindowedSinc:
  458.     FilterSupport = FilterRadius;
  459.     SetFilterSupport( FilterSupport );
  460.     window = 1;
  461.     BuildCdf( WindowedSincFunction );
  462.     break;
  463.   case Cubic:
  464.     FilterSupport = 2.0;
  465.     CubicFilterSetup( B, C );
  466.     break;
  467.   default:
  468.     fprintf( stderr, "Filter type: %d not available\n", FilterFunction );
  469.     exit(1);
  470.   }
  471.  
  472.   /* Start with the root node on the stack */
  473.   nStack = 0;
  474.   NodeStack[nStack] = Root;
  475.   BoundaryStack[nStack] = RootBoundary;
  476.   levelStack[nStack++] = 0;
  477.  
  478.   while (nStack > 0) {
  479.  
  480.     /* Get the next node to process from the stack */
  481.     Node = NodeStack[--nStack];
  482.     Boundary = BoundaryStack[nStack];
  483.     level    = levelStack[nStack];
  484.     
  485.     /* Follow the tree down a leaf, pushing the side branches.
  486.      */
  487.     while (! IsLeaf(Node)) {
  488.  
  489.       SplitBoundary( TRUE, Node->splitDirection, Node->splitValue,
  490.             &Boundary, &LeftBoundary );
  491.       if (Node->Left) {
  492.     NodeStack[nStack] = Node->Left;
  493.     BoundaryStack[nStack] = LeftBoundary;
  494.     levelStack[nStack++] = level+1;
  495.       }
  496.       SplitBoundary( FALSE, Node->splitDirection, Node->splitValue, 
  497.             &Boundary, &RightBoundary );
  498.       if (Node->Right) {
  499.     if (nStack >= MAX_DEPTH) {
  500.       fprintf( stderr, "Filter stack overflow!\n"  );
  501.       abort();
  502.     }
  503.     NodeStack[nStack] = Node->Right;
  504.     BoundaryStack[nStack] = RightBoundary;
  505.     levelStack[nStack++] = level+1;
  506.       }
  507.  
  508.       /* Get the next node to process from the stack */
  509.       Node = NodeStack[--nStack];
  510.       Boundary = BoundaryStack[nStack];
  511.       level    = levelStack[nStack];
  512.     }
  513.   
  514.     if (Node == NULL)
  515.       continue;        /* Dead end, go back to the stack loop */
  516.   
  517.     if (Animating)
  518.       AnimateCheckInput();
  519.  
  520.     /* Map the sample into the screen coordinate system */
  521.     Xs = Xsize*(Node->Sample->Xs+1)/2.0 - 0.5;
  522.     Ys = Ysize*(Node->Sample->Ys+1)/2.0 - 0.5;
  523.     Boundary.xH = Xsize*(Boundary.xH+1)/2.0 - 0.5;
  524.     Boundary.xL = Xsize*(Boundary.xL+1)/2.0 - 0.5;
  525.     Boundary.yL = Ysize*(Boundary.yL+1)/2.0 - 0.5;
  526.     Boundary.yH = Ysize*(Boundary.yH+1)/2.0 - 0.5;
  527.  
  528.     /* Find the cell that the sample actually lies in */
  529.     c = (int) floor( Xs );
  530.     r = (int) floor( Ys );
  531.     if (r >= y_start &&  r <= y_stop  &&
  532.     c >= x_start &&  c <= x_stop )  
  533.       image_data[r-y_start][c-x_start].samplesInPixel++;
  534.  
  535.     /* Find the rectangle pixels affected by the sample cell */
  536.     rmn = (int) floor ( Boundary.yL - FilterRadius );
  537.     rmx = (int) ceil  ( Boundary.yH + FilterRadius );
  538.     cmn = (int) floor ( Boundary.xL - FilterRadius );
  539.     cmx = (int) ceil  ( Boundary.xH + FilterRadius );
  540.  
  541.     /* loop over all affected pixels */
  542.     for(r=rmn; r<rmx; r++) {
  543.       if (r < y_start || r > y_stop) continue;
  544.       for(c=cmn; c<cmx; c++) {
  545.     if (c < x_start || c > x_stop) continue;
  546.  
  547.     /* Translate coordinates so that this affected pixel is at (0,0) and
  548.        so that the filter radius is FilterSupport units long 
  549.      */
  550.     yL = (FilterSupport * (Boundary.yL - r)) / FilterRadius;
  551.     yH = (FilterSupport * (Boundary.yH - r)) / FilterRadius;
  552.     xL = (FilterSupport * (Boundary.xL - c)) / FilterRadius;
  553.     xH = (FilterSupport * (Boundary.xH - c)) / FilterRadius;
  554.  
  555.     /* Lookup the filter values */
  556.     if (FilterFunction == Cubic) {
  557.       Weight = IntegralLookup( xL, yL, xH, yH );
  558.     } else {
  559.       Weight = ( CdfLookup( xH ) - CdfLookup( xL) ) *
  560.                ( CdfLookup( yH ) - CdfLookup( yL) );
  561.     }
  562.     pixel = &(image_data[r-y_start][c-x_start]);
  563.     clr   = &(Node->Sample->Ave);
  564.     pixel->weight += Weight;
  565.     pixel->clr.r += Weight*clr->r;
  566.     pixel->clr.g += Weight*clr->g;
  567.     pixel->clr.b += Weight*clr->b;
  568.     pixel->clr.a += Weight*clr->a;
  569.       }
  570.     }
  571.   }
  572. }
  573.  
  574.  
  575.  
  576.  
  577. ReSampleImage( Root, radius, stddevs )
  578. HierarchicalRegion *Root;
  579. double radius;
  580. double stddevs;
  581. {
  582.   int xpixel, ypixel;        /* The coordinate of the current pixel */
  583.   FrameBufferType *FrameBuffer = NULL;  /*WFF frame buffer for image*/
  584.   FrameBufferType *NumRaysFB = NULL;    /* for number-of-rays image */
  585.   FloatColor pixelclr;
  586.   static char *RGBABands = "RGBA";
  587.   static char *IBandOnly = "I";
  588.   char Value[ValueLength];
  589.   int x, y;
  590.   short inpixelcount;
  591.   float totalW;
  592.   float Radius;
  593.  
  594.   extern int y_stop, y_start, x_stop, x_start;
  595.   extern double Aspect;
  596.   extern FILE *image, *NumRaysImage;
  597.   extern int Xsize, Ysize;
  598.   int maxinpixel = 0;
  599.   int Dims[2];
  600.   Pixel *pixel;
  601.  
  602.   /* Allocate memory for the image */
  603.   Dims[0] = y_stop - y_start + 1;
  604.   Dims[1] = x_stop - x_start + 1;
  605.   image_data = (Pixel **) mallocNd( 2, Dims, sizeof(Pixel) );
  606.   assert( image_data );
  607.   bzero( image_data[0], Dims[0]*Dims[1]*sizeof(Pixel) );
  608.  
  609.   /* Traverse the tree and filter the samples */
  610.   Filter( Root, RootBoundary );
  611.  
  612.  
  613.   /* Set up the frame buffer */
  614.  
  615.   OpenFB(&FrameBuffer);
  616.   SetBounds(FrameBuffer, y_start, x_start, y_stop, x_stop);
  617.   SetColorSystem(FrameBuffer, RGBABands, BitsPerBand);
  618.   sprintf(Value, "%f", Aspect);
  619.   SetDescriptor(FrameBuffer, "AspectRatio", Value);
  620.   strcpy( Value, "RLE" );
  621.   SetDescriptor(FrameBuffer, "Encoding", Value );
  622.  
  623.   PassImageOut(image, FrameBuffer);
  624.  
  625.   /*------------------------------------------------------------*/
  626.   /* Initialize the number of rays image if it was requested.   */
  627.   /*------------------------------------------------------------*/
  628.   if (NumRaysImage != NULL) {
  629.     OpenFB(&NumRaysFB);
  630.     SetBounds(NumRaysFB, y_start, x_start, y_stop, x_stop);
  631.     SetColorSystem(NumRaysFB, IBandOnly, 8);
  632.     sprintf(Value, "%f", Aspect);
  633.     SetDescriptor(NumRaysFB, "AspectRatio", Value);
  634.     strcpy( Value, "RLE");
  635.     SetDescriptor( NumRaysFB, "Encoding", Value );
  636.     PassImageOut(NumRaysImage, NumRaysFB);
  637.   }
  638.  
  639.   fprintf(stderr,
  640. "Xsize: %d, Ysize: %d, x_start: %d, x_stop: %d, y_start: %d, y_stop: %d\n",
  641.       Xsize, Ysize, x_start, x_stop, y_start, y_stop);
  642.  
  643.   maxinpixel = 0;
  644.   for(y=y_start; y<=y_stop; y++) {
  645.     for(x=x_start; x<=x_stop; x++) {
  646.       pixel = &(image_data[y-y_start][x-x_start]);
  647.       if (pixel->weight != 0) {
  648.     pixel->clr.r /= pixel->weight;
  649.     pixel->clr.g /= pixel->weight;
  650.     pixel->clr.b /= pixel->weight;
  651.       }
  652.       Output( FrameBuffer, x, y,  pixel->clr );
  653.       if (NumRaysImage) 
  654.     NextPixelOut(NumRaysFB, &pixel->samplesInPixel);
  655.       if (pixel->samplesInPixel > maxinpixel)
  656.     maxinpixel = pixel->samplesInPixel;
  657.     }
  658.   }
  659.  
  660.   CloseFB(&FrameBuffer);
  661.   if (NumRaysImage)
  662.     CloseFB(&NumRaysFB);
  663.   fprintf( stderr, "Maximum samples in pixel: %d\n", maxinpixel );
  664.  
  665. }
  666.  
  667. void ParseCommandLine(argc, argv)
  668. int argc;
  669. char **argv;
  670. {
  671.   int nrays = 0;
  672.  
  673.   /* Jump over the command name */
  674.   argc--; argv++;
  675.  
  676.   /* Parse the arguments */
  677.   while (argc > 0) {
  678.     if (argv[0][0] != '-') {
  679.       fprintf(stderr, "Error in command line.\n");
  680.       Usage();
  681.       exit(-1);
  682.     }
  683.  
  684.     /* Determine which switch, if any. */
  685.  
  686.     /**** Please keep these lexographically sorted! *****/
  687.     switch (argv[0][1]) {
  688.     case 'A':
  689.       /* Set the aspect ratio of pixels (width/height) */
  690.       Aspect = (double) atof(argv[1]);
  691.       argc -= 2; argv += 2;
  692.       break;
  693.     case 'B':
  694.       Animating = TRUE;
  695.       argc--; argv++;
  696.       break;
  697.     case 'a':
  698.       /* Set anitaliaing parameters */
  699.       switch (argv[0][2]) {
  700.       case 'B':
  701.     B = atof(argv[1]);
  702.     argc--; argv++;
  703.     break;
  704.       case 'C':
  705.     CutOff = atof(argv[1]);    /* For Sinc filter option */
  706.                                 /* cutoff frequency in cycles-per-pixel */
  707.     C = CutOff;        /* For Cubic filter option */
  708.     argc--; argv++;
  709.     break;
  710.       case 'F': {
  711.     /* Set the filter function */
  712.     char *string = argv[1];
  713.     if (strcmp(string,"Gaussian") == 0) 
  714.       FilterFunction = Gaussian;
  715.     else if (strcmp(string, "Sinc") == 0) {
  716.       FilterFunction = Sinc;
  717.       if (CutOff < 0.0) {
  718.         CutOff = 1/2;  /* Nyquist frequency: half of sample rate */
  719.       }
  720.     } else if (strcmp(string, "WindowedSinc") == 0) {
  721.       FilterFunction = WindowedSinc;
  722.       if (CutOff < 0.0) {
  723.         CutOff = 1/2;  /* Nyquist frequency: half of sample rate */
  724.       }
  725.     } else if (strcmp(string, "Cubic") == 0)
  726.       FilterFunction = Cubic;
  727.     else {
  728.       fprintf( stderr, "Invalid filter type: %s\n\n", string );
  729.       Usage();
  730.       exit(1);
  731.     }
  732.     argc -= 1; argv += 1;
  733.     break;
  734.       }
  735.       case 'f':
  736.     /* Set the size of the support of the filter (in pixels) */
  737.     FilterRadius = (double) atof(argv[1]);
  738.     argc -= 1; argv += 1;
  739.     break;
  740.       case 'n':
  741.     {
  742.       char *nraysname = argv[1];
  743.  
  744.       /* Generate the nrays.im image */
  745.       if ((NumRaysImage = fopen( nraysname, "w")) == NULL) {
  746.         fprintf(stderr, "Can't open nrays.im image file.\n");
  747.         Usage();
  748.         exit(1);
  749.       }
  750.       argc -= 1; argv += 1;
  751.       break;
  752.     }
  753.       case 's':
  754.     /* Split nodes adaptively */
  755.     AdaptiveSplits = TRUE;
  756.     break;
  757.       case 't':
  758.     /* Set the width (in stddevs) of the truncated Gaussian */
  759.     FilterSupport = (double) atof(argv[1]);
  760.     argc -=1; argv += 1;
  761.     break;
  762.       default:
  763.     fprintf(stderr, "Unknown antialiasing option: %c\n", argv[0][2]);
  764.     Usage();
  765.     exit(1);
  766.       }
  767.       argc -= 1; argv += 1;
  768.       break;
  769.     case 'b':
  770.       BitsPerBand = atoi( argv[1] );
  771.       argc -= 2; argv +=2;
  772.       MaxOutBandValue = (1 << BitsPerBand) - 1;
  773.       break;
  774.     case 'i':
  775.       /* Set the image file name */
  776.       filename = argv[1];
  777.       argc -= 2; argv += 2;
  778.       break;
  779.     case 'm':
  780.       /* Set the maximum value expect in a color band */
  781.       MaxInBandValue = atof( argv[1] );
  782.       argc -= 2; argv += 2;
  783.       break;
  784.     case 's':
  785.       /* Set the start and stop scan lines    */
  786.       y_start = atoi(argv[1]);
  787.       y_stop  = atoi(argv[2]);
  788.       argc -= 3; argv += 3;
  789.       break;
  790.     case 'x':
  791.       /* Set the x resolution */
  792.       Xsize = atoi(argv[1]);
  793.       x_start = 0;
  794.       x_stop = Xsize-1;
  795.       argc -=2; argv += 2;
  796.       break;
  797.     case 'y':
  798.       /* Set the y resolution */
  799.       Ysize = atoi(argv[1]);
  800.       y_start = 0;
  801.       y_stop = Ysize-1;
  802.       argc -=2; argv += 2;
  803.       break;
  804.     case 'w':
  805.       /* Set the entire scanning window */
  806.       x_start = atoi(argv[1]);
  807.       x_stop  = atoi(argv[2]);
  808.       y_start = atoi(argv[3]);
  809.       y_stop  = atoi(argv[4]);
  810.       argc -= 5; argv += 5;
  811.       break;
  812.     default:
  813.       fprintf(stderr, "bad option: %s\n", argv[0]);
  814.       Usage();
  815.       exit(-1);
  816.     }
  817.   }
  818. }
  819.  
  820.  
  821.  
  822. main( argc, argv)
  823. int argc;
  824. char *argv[];
  825. {
  826.   RectType Boundary;        /* Boundary of interesting region */
  827.  
  828.   ParseCommandLine(argc, argv);
  829.  
  830.   if (filename == NULL)
  831.     image = stdout;
  832.   else
  833.     image = fopen( filename, "w" );
  834.   if (image == NULL) 
  835.     {
  836.       fprintf( stderr, "Unable to open output file\n" );
  837.       exit(1);
  838.     }
  839.  
  840.   PixelSizeX = 2.0/((double) Xsize);
  841.   PixelSizeY = 2.0/((double) Ysize);
  842.  
  843.   Boundary.xL = -1.0 + x_start*PixelSizeX - PixelSizeX*FilterRadius;
  844.   Boundary.xH = -1.0 + (x_stop+1)*PixelSizeX  + PixelSizeX*FilterRadius;
  845.   Boundary.yL = -1.0 + y_start*PixelSizeY - PixelSizeY*FilterRadius;
  846.   Boundary.yH = -1.0 + (y_stop+1)*PixelSizeY  + PixelSizeY*FilterRadius;
  847.  
  848.   if (Animating) {
  849.     Animate("");
  850.     AnimateInitialize( Boundary.xL, Boundary.yL, Boundary.xH, Boundary.yH );
  851.   }
  852.  
  853.  
  854.   /* Read the samples in */
  855.   ReadSampleFile( stdin, &Boundary );
  856.  
  857.   /* Filter and resample at the screen resolution */
  858.   ReSampleImage( &Root, FilterRadius, FilterSupport );
  859.  
  860.   if (Animating)
  861.     AnimateExit();
  862.  
  863.   return 0;
  864.  
  865. }
  866.  
  867.  
  868. Output(FrameBuffer, x, y, c)
  869. FrameBufferType *FrameBuffer;
  870. int x, y;
  871. FloatColor c;
  872. {
  873.   int red, green, blue, alpha;
  874.   unsigned short pixel[4];
  875.  
  876.   if ((red = (int) (c.r/MaxInBandValue * MaxOutBandValue  + 0.5)) 
  877.       > MaxOutBandValue)    red = MaxOutBandValue;
  878.   if (red < 0) red = 0;
  879.  
  880.   if ((green = (int) (c.g/MaxInBandValue * MaxOutBandValue + 0.5)) 
  881.       > MaxOutBandValue)  green = MaxOutBandValue;
  882.   if (green < 0) green = 0;
  883.  
  884.   if ((blue = (int) (c.b/MaxInBandValue * MaxOutBandValue + 0.5)) 
  885.       > MaxOutBandValue)   blue = MaxOutBandValue;
  886.   if (blue < 0) blue = 0;
  887.  
  888.   if ((alpha = (int) (c.a * MaxOutBandValue + 0.5)) 
  889.       >MaxOutBandValue)  alpha = MaxOutBandValue;
  890.   if (alpha < 0) alpha = 0;
  891.  
  892.   pixel[0] = red;
  893.   pixel[1] = green;
  894.   pixel[2] = blue;
  895.   pixel[3] = alpha;
  896.  
  897.   NextPixelOut(FrameBuffer, pixel);
  898.  
  899. }
  900.  
  901. Usage()
  902. {
  903.  
  904.   fprintf( stderr, "\
  905. Usage: Reconstruct [options].\n\
  906. Options are:\n\
  907. -aB val -aC val           \"Mitchell\" parameters for cubic filter.  \n\
  908. -aC cutoff                Cutoff wavelength for Sinc filter. In cycles/pixel.\n\
  909. -aF [ Cubic       |          Filter function \n\
  910.       Gaussian    |   \n\
  911.       Sinc        |   \n\
  912.       WindowedSinc   ] \n\
  913. -af radius                Filter radius (in pixels). \n\
  914. -an filename              Generate an image showing number of samples.\n\
  915. -as                       Split nodes adaptively.\n\
  916. -at support               Width of gaussian filter (in stddevs).\n\
  917. \n\
  918. -A aspect                 Set aspect ration\n\
  919. -i filename               Set result file name\n\
  920. -x xsize                  Set the resolution of the image in the x direction.\n\
  921. -y ysize                  Set the resolution of the image in the y direction.\n\
  922. -s y_start x_start        Set the start and stop scan lines\n\
  923. -m value                  Set the maximum value expected in a color band\n\
  924. -w xmin xmax ymin ymax    Set the entire scan window\n\
  925. ");
  926. }
  927.  
  928.